packages = c('rgdal', 'sf', 'tmap', 'tidyverse', 'sp', 'rgeos','maptools', 'raster', 'spatstat', 'tmaptools', 'spdep', 'OpenStreetMap', 'ggpubr', 'SpatialPosition', 'SpatialAcc', 'dplyr', 'shinycssloaders','plotly', 'shinythemes')
for (p in packages){
if (!require(p, character.only = T)){
install.packages(p)
}
library(p, character.only = T)
}
# Population Data
# columns = planning_area, subzone, elderly_count, total_count
popdata <- read_csv('data/planning-area-subzone-age-group-sex-and-type-of-dwelling-june-2011-2019.csv')
popdata2019 <- popdata %>%
filter(year == 2019) %>%
# filter(age_group == "65_to_69" | age_group == "70_to_74" | age_group == '75_to_79' | age_group == '80_to_84' | age_group == '85_to_89' | age_group == '90_and_over') %>%
group_by(planning_area, subzone, age_group) %>%
summarise(count = sum(resident_count)) %>%
ungroup() %>%
spread(age_group, count) %>%
mutate(elderly_count = `65_to_69` + `70_to_74` + `75_to_79` + `80_to_84` + `85_to_89` + `90_and_over`) %>%
mutate(total_count = rowSums(.[3:21])) %>%
dplyr::select(planning_area, subzone, elderly_count, total_count)
popdata2019 <- mutate_at(popdata2019, .vars = c("subzone", "planning_area"), .funs=toupper)
# Eldercare Data
eldercare_sf <- st_read(dsn='data', layer='ELDERCARE')
## Reading layer `ELDERCARE' from data source `D:\IS415Team2\codes\data' using driver `ESRI Shapefile'
## Simple feature collection with 133 features and 18 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 14481.92 ymin: 28218.43 xmax: 41665.14 ymax: 46804.9
## projected CRS: SVY21
eldercare_sf <- eldercare_sf %>%
mutate(label = "Eldercare centres")
eldercare_sf <- st_transform(eldercare_sf, 3414)
st_crs(eldercare_sf)
## Coordinate Reference System:
## User input: EPSG:3414
## wkt:
## PROJCRS["SVY21 / Singapore TM",
## BASEGEOGCRS["SVY21",
## DATUM["SVY21",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4757]],
## CONVERSION["Singapore Transverse Mercator",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",1.36666666666667,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",103.833333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",1,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",28001.642,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",38744.572,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["northing (N)",north,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["easting (E)",east,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["Singapore"],
## BBOX[1.13,103.59,1.47,104.07]],
## ID["EPSG",3414]]
eldercare_sp <- as_Spatial(eldercare_sf)
eldercare_spatialpoint <- as(eldercare_sp, "SpatialPoints")
# Silver Infocomm Data
infocomm_sf <- st_read(dsn='data', layer='SILVERINFOCOMM')
## Reading layer `SILVERINFOCOMM' from data source `D:\IS415Team2\codes\data' using driver `ESRI Shapefile'
## Simple feature collection with 101 features and 18 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 12024.87 ymin: 28385.78 xmax: 41300.53 ymax: 47763.05
## projected CRS: SVY21
infocomm_sf <- infocomm_sf %>%
mutate(label = "Silver Infocomm Junc")
infocomm_sf <- st_transform(infocomm_sf, 3414)
st_crs(infocomm_sf)
## Coordinate Reference System:
## User input: EPSG:3414
## wkt:
## PROJCRS["SVY21 / Singapore TM",
## BASEGEOGCRS["SVY21",
## DATUM["SVY21",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4757]],
## CONVERSION["Singapore Transverse Mercator",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",1.36666666666667,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",103.833333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",1,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",28001.642,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",38744.572,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["northing (N)",north,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["easting (E)",east,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["Singapore"],
## BBOX[1.13,103.59,1.47,104.07]],
## ID["EPSG",3414]]
infocomm_sp <- as_Spatial(infocomm_sf)
infocomm_spatialpoint <- as(infocomm_sp, "SpatialPoints")
# Chas Clinic Data
chas_sf <- st_read(dsn='data/chas-clinics-kml.kml')
## Reading layer `MOH_CHAS_CLINICS' from data source `D:\IS415Team2\codes\data\chas-clinics-kml.kml' using driver `KML'
## Simple feature collection with 1167 features and 2 fields
## geometry type: POINT
## dimension: XYZ
## bbox: xmin: 103.5818 ymin: 1.016264 xmax: 103.9903 ymax: 1.456037
## z_range: zmin: 0 zmax: 0
## geographic CRS: WGS 84
chas_sf <- chas_sf %>%
mutate(label = "Chas Clinics") %>%
mutate(capacity = 1)
chas_sf <- st_transform(chas_sf, 3414)
st_crs(chas_sf)
## Coordinate Reference System:
## User input: EPSG:3414
## wkt:
## PROJCRS["SVY21 / Singapore TM",
## BASEGEOGCRS["SVY21",
## DATUM["SVY21",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4757]],
## CONVERSION["Singapore Transverse Mercator",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",1.36666666666667,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",103.833333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",1,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",28001.642,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",38744.572,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["northing (N)",north,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["easting (E)",east,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["Singapore"],
## BBOX[1.13,103.59,1.47,104.07]],
## ID["EPSG",3414]]
chas_sp <- as_Spatial(chas_sf)
chas_spatialpoint <- as(chas_sp, "SpatialPoints")
# MP14 Subzone Data
mpsz_sf <- st_read(dsn='data', layer='MP14_SUBZONE_WEB_PL')
## Reading layer `MP14_SUBZONE_WEB_PL' from data source `D:\IS415Team2\codes\data' using driver `ESRI Shapefile'
## Simple feature collection with 323 features and 15 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
## projected CRS: SVY21
mpsz_sf <- mpsz_sf %>%
dplyr::select(SUBZONE_N, PLN_AREA_N, REGION_N, X_ADDR, Y_ADDR, SHAPE_Leng, SHAPE_Area, geometry)
mpsz_sf <- st_transform(mpsz_sf, 3414)
st_crs(mpsz_sf)
## Coordinate Reference System:
## User input: EPSG:3414
## wkt:
## PROJCRS["SVY21 / Singapore TM",
## BASEGEOGCRS["SVY21",
## DATUM["SVY21",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4757]],
## CONVERSION["Singapore Transverse Mercator",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",1.36666666666667,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",103.833333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",1,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",28001.642,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",38744.572,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["northing (N)",north,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["easting (E)",east,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["Singapore"],
## BBOX[1.13,103.59,1.47,104.07]],
## ID["EPSG",3414]]
eldercare_sf <- st_join(eldercare_sf, mpsz_sf, join=st_intersects)
infocomm_sf <- st_join(infocomm_sf, mpsz_sf, join=st_intersects)
chas_sf <- st_join(chas_sf, mpsz_sf, join=st_intersects)
# HDB Data
hdb <- read_csv('data/hdb_data.csv')
hdb <- hdb %>%
mutate(total_count=rowSums(.[2:11])) %>%
dplyr::select(Postcode, X, Y, Latitude, Longitude, total_count)
hdb$Postcode <- as.character(hdb$Postcode)
hdb_sf <- st_as_sf(hdb, coords=c('X', 'Y'), crs='EPSG:3414')
# SG Coastal Outline
sg <- readOGR(dsn = "data", layer="CostalOutline")
## OGR data source with driver: ESRI Shapefile
## Source: "D:\IS415Team2\codes\data", layer: "CostalOutline"
## with 60 features
## It has 4 fields
sg_spatialpoint <- as(sg, "SpatialPolygons")
# Compute elderly_density in mpsz
mpsz_sf <- left_join(mpsz_sf, popdata2019, by=c('PLN_AREA_N' = 'planning_area', 'SUBZONE_N' = 'subzone'))
mpsz_sf <- mpsz_sf %>%
dplyr::select(SUBZONE_N, PLN_AREA_N, REGION_N, elderly_count, total_count, geometry) %>%
mutate(elderly_proportion = elderly_count / total_count) %>%
mutate_if(is.numeric, ~replace(., is.nan(.), 0))
mpsz_sp <- as_Spatial(mpsz_sf)
mpsz_spatialpoint <- as(mpsz_sp, "SpatialPolygons")
# Add count of facilities per subzone into mpsz_demand
mpsz_demand <- mpsz_sf
mpsz_demand$`chas_count` <- lengths(st_intersects(mpsz_sf,chas_sf))
mpsz_demand$`eldercare_count` <- lengths(st_intersects(mpsz_sf,eldercare_sf))
mpsz_demand$`infocomm_count` <- lengths(st_intersects(mpsz_sf,infocomm_sf))
# Join HDB with MP14 Data
hdb_mpsz <- st_join(mpsz_sf, hdb_sf, join=st_intersects) %>%
dplyr::select(SUBZONE_N, PLN_AREA_N, REGION_N, elderly_proportion, Postcode, total_count.y, geometry) %>%
rename(resident_count = total_count.y) %>%
mutate(elderly_count = resident_count * elderly_proportion) %>%
mutate_if(is.numeric, ~replace(., is.na(.), 0))
#Finding area of each subzone
mpsz_demand$Area <- mpsz_demand %>%
st_area()
#Calculating density of Chas, Eldercare, Infocomm and Elderly Pop
mpsz_demand <- mpsz_demand %>%
mutate(`Elderly_Density` = (`elderly_count`/ Area) * 1000000) %>%
mutate(`Chas_Density` = `elderly_count` / `chas_count`) %>%
mutate(`Eldercare_Density` = `elderly_count` / `eldercare_count`) %>%
mutate(`Infocomm_Density` = `elderly_count` / `infocomm_count`)
#Converting mpsz_demand to SP object
mpsz_demand_sp <- as_Spatial(mpsz_demand)
mpsz_demand_spatialpoint <- as(mpsz_demand_sp, "SpatialPolygons")
#reading OSM
sg_osm <- read_osm(mpsz_spatialpoint, ext=1.3)
#Extracting unique planning areas
uniqPlanningAreas <- mpsz_demand[2]
uniqPlanningAreas <- st_set_geometry(uniqPlanningAreas, NULL)
uniqPlanningAreas <- unique(uniqPlanningAreas)
Quadrat Test
CHAS Clinics
owin <- as(sg, "owin")
spatialpoint <- as(sg, "SpatialPolygons")
chas_ppp <- as(chas_spatialpoint, "ppp")
chas_ppp_jit <- rjitter(chas_ppp, retry=TRUE, nsim=1, drop=TRUE)
ppp <- chas_ppp_jit[owin]
ppp.km <- rescale(ppp, 1000, "km")
qt <- quadrat.test(ppp.km, nx=20, ny=15)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
qt
##
## Chi-squared test of CSR using quadrat counts
##
## data: ppp.km
## X2 = 1976.7, df = 184, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 185 tiles (irregular windows)
plot(ppp.km)
plot(qt, add=TRUE, cex=.1)
quadrat.test(ppp.km, nx=20, ny=15, method='M', nsim=2999)
##
## Conditional Monte Carlo test of CSR using quadrat counts
## Test statistic: Pearson X2 statistic
##
## data: ppp.km
## X2 = 1976.7, p-value = 0.0006667
## alternative hypothesis: two.sided
##
## Quadrats: 185 tiles (irregular windows)
Eldercare
owin <- as(sg, "owin")
spatialpoint <- as(sg, "SpatialPolygons")
eldercare_ppp <- as(eldercare_spatialpoint, "ppp")
eldercare_ppp_jit <- rjitter(eldercare_ppp, retry=TRUE, nsim=1, drop=TRUE)
ppp <- eldercare_ppp_jit[owin]
ppp.km <- rescale(ppp, 1000, "km")
qt <- quadrat.test(ppp.km, nx=20, ny=15)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
qt
##
## Chi-squared test of CSR using quadrat counts
##
## data: ppp.km
## X2 = 495.78, df = 184, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 185 tiles (irregular windows)
plot(ppp.km)
plot(qt, add=TRUE, cex=.1)
quadrat.test(ppp.km, nx=20, ny=15, method='M', nsim=2999)
##
## Conditional Monte Carlo test of CSR using quadrat counts
## Test statistic: Pearson X2 statistic
##
## data: ppp.km
## X2 = 495.78, p-value = 0.007333
## alternative hypothesis: two.sided
##
## Quadrats: 185 tiles (irregular windows)
Silver Infocomm Junctions
owin <- as(sg, "owin")
spatialpoint <- as(sg, "SpatialPolygons")
infocomm_ppp <- as(infocomm_spatialpoint, "ppp")
infocomm_ppp_jit <- rjitter(infocomm_ppp, retry=TRUE, nsim=1, drop=TRUE)
ppp <- infocomm_ppp_jit[owin]
ppp.km <- rescale(ppp, 1000, "km")
qt <- quadrat.test(ppp.km, nx=20, ny=15)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
qt
##
## Chi-squared test of CSR using quadrat counts
##
## data: ppp.km
## X2 = 317.48, df = 184, p-value = 7.149e-09
## alternative hypothesis: two.sided
##
## Quadrats: 185 tiles (irregular windows)
plot(ppp.km)
plot(qt, add=TRUE, cex=.1)
quadrat.test(ppp.km, nx=20, ny=15, method='M', nsim=2999)
##
## Conditional Monte Carlo test of CSR using quadrat counts
## Test statistic: Pearson X2 statistic
##
## data: ppp.km
## X2 = 317.48, p-value = 0.05133
## alternative hypothesis: two.sided
##
## Quadrats: 185 tiles (irregular windows)
Gi* Statistics - Hotspot & Coldspot Analysis
CHAS
#coverting to sp polygons
mpsz_demand_selected_sp <- as_Spatial(mpsz_demand)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum SVY21 in CRS definition
coords <- coordinates(mpsz_demand_selected_sp)
k1 <- knn2nb(knearneigh(coords))
k1dists <- unlist(nbdists(k1, coords, longlat=FALSE))
maxdist <- ceiling(max(k1dists))
#Calculating d fixed-distance weight matrix
dnb <- dnearneigh(coordinates(mpsz_demand_selected_sp), 0, maxdist, longlat = FALSE)
#plot(mpsz_demand_selected_sp, border = 'lightgrey')
#plot(dnb, coordinates(mpsz_demand_selected_sp), add=TRUE, col='red')
dnb_lw <- nb2listw(dnb, style = 'B')
fips <- order(mpsz_demand_selected_sp$SUBZONE_N)
gi.fixed <- localG(mpsz_demand_selected_sp$chas_count, dnb_lw)
gi.fixed
## [1] -1.097895747 -1.750276041 -1.224729244 -1.649028543 -2.064212720
## [6] -2.104086969 -1.901554875 -1.369558330 -2.593432159 -2.101292745
## [11] -2.747198138 -1.546764134 -0.443062227 -1.503022100 -1.619944188
## [16] -2.883406367 -1.032641641 -0.729049655 -1.034345429 -2.321395295
## [21] -1.712554693 -1.380816216 -1.377488708 -2.354621229 -1.945165127
## [26] -1.391393537 -1.480648481 -1.893688362 -2.471743794 -2.791361849
## [31] -1.341978208 -1.900254240 -1.853839054 -1.577809366 -2.066549368
## [36] -1.635487565 -1.350746913 -1.433679257 -1.406179580 -1.903140385
## [41] -0.954469665 -1.845231592 -1.176629161 0.597843958 -0.951053930
## [46] -0.578373599 -1.238172440 -1.201281267 -2.299024629 2.470425585
## [51] -1.925531529 -1.102744323 -1.008256238 -0.222250411 -1.232462559
## [56] -2.129614118 -0.956845194 -0.109051056 -0.061443232 -1.012196402
## [61] -2.099018138 -1.401456246 -1.924088906 -1.804893203 -0.131215798
## [66] 0.208445774 -1.785318927 0.607510237 1.244966241 -1.766319485
## [71] -2.214459555 -2.142616072 -0.773761038 -1.198891758 -1.547837882
## [76] -0.519272574 2.963033961 -1.518153429 -1.708842703 -0.884783701
## [81] -0.524114540 1.206227662 -1.004389124 -0.460550974 0.170551313
## [86] -1.448377090 -1.275613027 -1.453325491 -2.510166044 -0.775924961
## [91] -0.781738019 3.867165618 -1.243749916 -1.364920993 -1.561403798
## [96] -0.164861841 -0.400156963 -0.663477304 0.091373939 -0.181244996
## [101] -2.042811789 -1.096868369 2.263711543 -0.610344516 1.035179183
## [106] -0.461078963 -0.835681162 3.359510264 0.403639797 -0.092678780
## [111] -1.760848295 -1.149713099 1.929301293 -0.691916596 -0.689440003
## [116] -0.991327362 3.509434389 2.492025187 -1.783106419 -0.020486411
## [121] -0.316626207 -0.748973304 -0.640682356 -1.340390173 -1.759022989
## [126] 0.298219769 -0.487173476 -1.113412729 -1.707779143 -1.560725383
## [131] 0.150243991 1.006615778 0.444707511 -1.215477121 -0.700423409
## [136] -0.753898934 3.692363911 1.656606566 0.519779890 -1.767771086
## [141] 0.581051086 0.331250221 -0.744045789 -0.854269586 0.727213838
## [146] -0.752342166 -0.353957904 -0.135514849 4.816891117 -0.484375283
## [151] -0.665464144 -0.660591092 3.650728407 -0.543479673 -0.823613899
## [156] -0.804595171 2.226565860 -0.433703572 0.341962401 -0.605707770
## [161] 0.351531512 1.996550070 1.464655849 -1.770583299 -0.660028257
## [166] -0.615901112 2.701726336 -1.030122195 -1.585272657 2.726082529
## [171] -0.997068570 -0.378124714 -0.667194145 -0.656505414 -0.631230573
## [176] 2.116287814 -0.779633114 -1.191791490 4.438419343 1.797010574
## [181] 0.189959224 3.816613245 -0.553424696 4.109907249 -0.377496888
## [186] 0.410182962 -0.446047300 -0.951574046 2.830403238 3.762186860
## [191] -1.539150633 -0.407622050 -1.033057370 0.611124624 -1.118413741
## [196] -1.049475304 2.104382260 0.542541079 2.915818858 -0.024313786
## [201] 3.307348449 5.154034976 5.059951346 4.959442232 1.848243676
## [206] -0.485231877 0.018343884 4.603710056 3.854585776 0.180737056
## [211] -0.970553433 -0.304491459 3.617087201 3.115111856 0.055083410
## [216] 3.646815817 1.737154398 2.964538659 0.259977405 0.229945954
## [221] 2.627725713 -0.209747847 -0.094832404 -0.175113724 2.041998578
## [226] 3.288052444 2.060931592 2.009458242 2.675943006 2.038377353
## [231] 3.558856844 2.845958537 -0.777569553 -0.371987115 3.831941212
## [236] 0.262704797 -0.438097967 1.415084478 2.763624792 1.068759475
## [241] 1.004496532 1.579760687 4.983842000 2.442161253 -0.482342250
## [246] 2.785630534 3.927556462 1.912012549 1.940969274 1.686096189
## [251] 2.450373117 -0.270167897 3.059914933 3.734143532 -0.165299297
## [256] 4.341412597 -0.561380812 3.547369056 0.218680739 3.427092814
## [261] 1.260432448 1.425562372 3.149260431 0.455004649 3.398295195
## [266] 1.651349191 1.848801277 -0.016203023 3.866855809 2.143393269
## [271] 1.812849398 2.886512051 1.697836081 3.812651631 1.120558764
## [276] 2.175207367 0.286056637 0.068167348 0.915814740 1.666586860
## [281] -0.450310423 -0.370101867 -0.277393718 1.166187367 0.619002233
## [286] 1.351780797 -0.420218122 -0.168659987 -0.257462928 -1.045887948
## [291] -0.670982576 -0.049261540 -0.383978403 -0.402871897 0.163324809
## [296] -0.578440303 -0.517922395 -0.151407320 -0.402871897 -0.214008828
## [301] -0.130580152 -0.729049655 -0.496288425 1.284387953 0.121629718
## [306] -0.440804399 -0.169385777 -0.565079733 1.045040426 1.324366000
## [311] 0.276377584 1.482616983 1.164900538 1.141327426 0.306010504
## [316] 1.334154568 0.642349641 -0.292232586 1.421088188 -0.106483069
## [321] -0.454036964 0.009388915 -0.155004360
## attr(,"gstari")
## [1] FALSE
## attr(,"call")
## localG(x = mpsz_demand_selected_sp$chas_count, listw = dnb_lw)
## attr(,"class")
## [1] "localG"
chas.gi <- cbind(mpsz_demand_selected_sp, as.matrix(gi.fixed))
names(chas.gi)[15] <- "gstat"
#colnames(chas.gi)[23] <- "gstat"
tm_shape(chas.gi) +
tm_fill(col = "gstat",
style = "pretty",
palette="-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## Variable(s) "gstat" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
Silver Infocomm
#coverting to sp polygons
mpsz_demand_selected_sp <- as_Spatial(mpsz_demand)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum SVY21 in CRS definition
coords <- coordinates(mpsz_demand_selected_sp)
k1 <- knn2nb(knearneigh(coords))
k1dists <- unlist(nbdists(k1, coords, longlat=FALSE))
maxdist <- ceiling(max(k1dists))
#Calculating d fixed-distance weight matrix
dnb <- dnearneigh(coordinates(mpsz_demand_selected_sp), 0, maxdist, longlat = FALSE)
#plot(mpsz_demand_selected_sp, border = 'lightgrey')
#plot(dnb, coordinates(mpsz_demand_selected_sp), add=TRUE, col='red')
dnb_lw <- nb2listw(dnb, style = 'B')
fips <- order(mpsz_demand_selected_sp$SUBZONE_N)
gi.fixed <- localG(mpsz_demand_selected_sp$infocomm_count, dnb_lw)
gi.fixed
## [1] -0.486612476 -1.376364058 -1.008135286 -1.205008596 -1.917588180
## [6] -1.835112563 -1.515746348 -1.179188873 -1.760654330 -1.924128231
## [11] -1.861020247 -1.234437804 0.149127773 -1.319506248 -1.319506248
## [16] -1.645565831 -0.551036959 -0.389034578 -0.882050612 -1.641850457
## [21] -1.064487068 -1.248137965 -1.467860028 -2.092487861 -1.771531599
## [26] -1.275772612 -1.318327319 -1.518589515 -2.288893189 -1.882044464
## [31] -1.290805083 -1.560752690 -1.738600475 -1.162327170 -1.936417346
## [36] -1.403206798 -1.063604526 -0.977752383 -1.150709925 -1.261388659
## [41] -0.934697775 -1.587150494 -1.063604526 0.421790648 -0.236677231
## [46] -0.953129260 -1.252335640 -1.008135286 -1.928944488 2.088328803
## [51] -1.374642034 -1.294917439 -1.294917439 0.003901357 -1.264393995
## [56] -1.524727012 -1.294917439 -0.076361780 -0.039843614 -1.294917439
## [61] -1.750221874 -1.190404776 -1.434030182 -1.487816736 -0.951978511
## [66] -0.013825444 -1.502838754 0.343619428 0.709842797 -1.319506248
## [71] -1.291053783 -1.768591752 -1.326050711 -0.114128677 -0.875378934
## [76] -1.176441705 3.233741962 -0.821940607 0.210542781 0.282507712
## [81] -1.261983618 0.239061681 -1.272460917 -1.133662674 0.063885256
## [86] 0.217861017 -1.136513876 -0.994864430 -1.486036168 -1.443230893
## [91] 1.153093188 3.547706765 -0.049266883 0.132227737 -1.319506248
## [96] -0.179296261 -0.973171863 -1.198144937 -0.127106267 -0.004347985
## [101] -1.684449627 -1.283432859 1.846875306 1.247657145 -0.090438037
## [106] 1.153093188 -1.604659789 2.882880156 -0.294967230 -0.085358807
## [111] -0.494449761 0.085349263 0.721540310 0.783652921 0.687368524
## [116] -0.793907776 3.211302257 2.115279564 -1.777480801 -0.143145261
## [121] -1.176441705 0.503994319 0.245843706 0.017271964 -1.292241092
## [126] 0.115521520 -1.390330359 0.838833982 -0.844843424 0.826178498
## [131] -0.080845811 0.168626291 -0.110211408 -0.465143267 0.877980656
## [136] -0.342860900 2.567701077 0.582915024 -0.145219920 0.115445788
## [141] -0.136905521 -0.180199522 0.958002207 0.658991860 -0.387461980
## [146] 0.828357832 1.286702686 -1.401314417 3.481010487 1.286702686
## [151] 0.997662419 0.795022020 3.179633314 1.048620418 -1.500787116
## [156] -0.819098952 1.215230212 1.095332948 -0.456505456 1.001351757
## [161] -0.223477078 0.942163267 0.304532611 -1.773714843 -1.411987426
## [166] 0.750348145 2.951299975 0.635832044 -1.005233344 1.723438657
## [171] 0.769360026 0.854653912 1.057173571 1.044125256 0.747455913
## [176] 1.025748504 1.000442492 -0.132419064 2.921121419 0.894749185
## [181] -0.369523103 2.726888733 1.011143727 2.812836350 1.170642824
## [186] -0.521022855 1.420865341 0.167252189 2.620998241 2.328978335
## [191] 0.300711918 1.187411662 0.861530148 0.012105359 -1.647632551
## [196] 0.583073371 1.120556466 -0.223477078 3.284272272 -1.167513969
## [201] 2.246487378 3.135308815 3.208296650 3.094768036 -0.993580210
## [206] 0.601204966 1.222362162 2.786944939 2.848951994 0.271836140
## [211] 0.371328639 -0.312680521 3.727028771 3.552917404 1.515209891
## [216] 2.966972941 -0.574739535 1.900328919 1.347524120 1.501708413
## [221] 1.733797638 0.050968122 0.857044104 0.857044104 -0.528908873
## [226] 2.425195734 -0.589227639 -0.449999441 2.276903291 -0.729821079
## [231] 2.533257431 3.331309421 0.893907436 -0.875378934 2.824346690
## [236] -0.867118969 0.866905546 -0.607031413 1.862051396 -0.739070529
## [241] -1.162496176 0.336873331 2.992349884 1.564633673 0.569271437
## [246] 3.546558700 2.763942285 -0.183343531 -0.473621867 -0.782606833
## [251] 0.102466886 1.321905465 1.936338196 3.057574722 1.297014958
## [256] 2.924314950 0.085349263 3.455600764 1.673941810 1.951219237
## [261] -0.977351573 -1.067251071 2.180102726 1.234348476 2.763942285
## [266] 0.509209071 1.158429208 0.565629412 2.809967322 -0.821940607
## [271] 1.512947431 1.943021784 -0.422786180 2.563991079 -1.487632051
## [276] -0.348196319 -0.944508035 -1.874206604 -1.621224918 -0.778358810
## [281] -0.145078137 -1.046703609 -0.672850225 -1.303424879 2.320956289
## [286] -0.828211667 -1.696496133 -0.951635669 0.683429518 -1.635542723
## [291] -1.854800629 -1.420266926 -1.768359077 -1.480074459 -1.133527045
## [296] -1.529900947 0.300711918 -1.754732685 -1.480074459 -1.377921639
## [301] -1.788188205 -0.389034578 -1.650802421 -0.982254895 -1.558554170
## [306] -0.926921824 -0.281736893 -1.800177349 -1.132824855 0.209159055
## [311] -1.090483138 1.066992242 0.371267883 -0.807216324 0.915479515
## [316] -0.642394728 -1.614114182 -0.573319219 1.273907405 -1.318381795
## [321] -1.529900947 -1.325485892 -1.377921639
## attr(,"gstari")
## [1] FALSE
## attr(,"call")
## localG(x = mpsz_demand_selected_sp$infocomm_count, listw = dnb_lw)
## attr(,"class")
## [1] "localG"
infocomm.gi <- cbind(mpsz_demand_selected_sp, as.matrix(gi.fixed))
names(infocomm.gi)[15] <- "gstat"
#colnames(infocomm.gi)[23] <- "gstat"
tm_shape(infocomm.gi) +
tm_fill(col = "gstat",
style = "pretty",
palette="-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## Variable(s) "gstat" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
Eldercare
#coverting to sp polygons
mpsz_demand_selected_sp <- as_Spatial(mpsz_demand)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum SVY21 in CRS definition
coords <- coordinates(mpsz_demand_selected_sp)
k1 <- knn2nb(knearneigh(coords))
k1dists <- unlist(nbdists(k1, coords, longlat=FALSE))
maxdist <- ceiling(max(k1dists))
#Calculating d fixed-distance weight matrix
dnb <- dnearneigh(coordinates(mpsz_demand_selected_sp), 0, maxdist, longlat = FALSE)
#plot(mpsz_demand_selected_sp, border = 'lightgrey')
#plot(dnb, coordinates(mpsz_demand_selected_sp), add=TRUE, col='red')
dnb_lw <- nb2listw(dnb, style = 'B')
fips <- order(mpsz_demand_selected_sp$SUBZONE_N)
gi.fixed <- localG(mpsz_demand_selected_sp$eldercare_count, dnb_lw)
gi.fixed
## [1] 2.093809425 1.318729340 1.698574954 0.614115327 0.674046348
## [6] 0.531687060 1.339159707 1.612148831 -0.444692848 -0.234421411
## [11] -0.971363307 1.115574781 2.067513040 1.372926746 1.372926746
## [16] -1.847163555 -0.618544314 -0.436695075 1.236712454 0.350096505
## [21] 0.653143550 1.917090489 2.067517888 -0.097971533 0.616302086
## [26] 1.824586859 1.771542473 0.950192536 -0.236360393 -0.823440130
## [31] 1.863194674 1.581660678 0.306931143 1.613853860 0.609955954
## [36] 1.523766969 1.682282507 1.787074608 1.682282507 1.592337298
## [41] 1.839819811 0.936735466 1.682282507 2.271290458 1.663511226
## [46] 2.194658797 1.833906951 1.818781797 -0.004908287 1.987394780
## [51] 0.559135889 2.139730569 2.139730569 2.547811491 1.868006495
## [56] 0.244764736 2.139730569 2.190859221 2.429170865 2.139730569
## [61] 0.699965243 1.931448968 1.130656597 0.356371372 2.241820468
## [66] 2.463870810 -0.702763791 2.324121450 1.289110801 0.465214467
## [71] 0.416195455 0.932081276 2.173898341 1.097832957 -0.982621316
## [76] 2.141172103 2.238246073 0.707037179 -1.027893317 -0.877331911
## [81] 2.159213134 2.615226252 2.309735762 2.191419676 2.529346174
## [86] -1.023010008 1.663511226 1.363207076 -1.668089965 2.109325218
## [91] -0.387470456 1.608045340 -0.921127080 -0.881287655 1.612148831
## [96] 2.377478410 2.309843311 2.140707983 2.275338852 1.991886291
## [101] 0.584862277 2.224534580 1.801634120 -0.766656152 2.228337860
## [106] -0.387470456 2.116696707 1.048958800 2.301445676 2.174039173
## [111] -1.903839227 -0.789987361 2.563745958 -0.822752724 -0.737394625
## [116] 1.469655842 2.066337705 1.474232292 0.106271036 2.377652095
## [121] 2.377652095 -0.893708693 -0.735288835 -0.856209401 -0.560933211
## [126] 2.245021117 2.363732232 -0.469815627 0.289407435 0.462478944
## [131] 2.301445676 2.343107686 2.209240004 0.618884550 -1.023957135
## [136] 1.113501705 1.545293926 2.145310185 1.892011266 -1.828348667
## [141] 2.036154379 2.242716270 -0.797139415 -0.491527619 1.854734050
## [146] -1.355266687 -0.897644064 1.998960252 2.691078771 -0.897644064
## [151] -0.769768776 -0.424766489 1.533931833 -0.604199052 1.803324827
## [156] 1.243647025 1.545360967 -0.571478882 1.696651548 -1.023957135
## [161] 2.352169617 1.631669049 1.768207936 1.145722142 2.028878949
## [166] -0.677763415 1.940268440 0.523336583 0.564511093 1.367536388
## [171] -0.692841005 -0.845748361 -0.971187610 -1.113864123 -1.138058248
## [176] 0.932081276 -1.095495196 -1.668089965 2.818601813 1.032592479
## [181] 1.852829350 0.734585952 -0.985950335 2.662697137 -0.864375684
## [186] 0.963594686 -0.581080557 -0.231627044 -0.421498973 2.055646099
## [191] -0.952934588 -1.180051357 -0.735288835 0.851317257 1.714492812
## [196] 0.076187155 1.425944264 2.239237880 0.864663507 2.072479360
## [201] 0.358418356 2.797607307 2.806962143 2.723733544 0.289407435
## [206] 0.071683671 -0.509826981 3.007235385 1.802136107 0.742773366
## [211] -0.874602082 1.613853860 0.717103125 0.560179418 -0.332568742
## [216] 0.266997709 0.693068610 2.083517689 0.089514422 0.378612796
## [221] 0.272253000 0.399168079 0.270875752 0.270875752 1.057629860
## [226] 2.374192653 1.302198555 1.083785797 0.325579368 -0.014506531
## [231] 2.532334922 0.316796962 -0.832014308 -0.982621316 1.755142250
## [236] 0.784945513 -0.877331911 0.968921740 2.025745544 0.776694173
## [241] 0.319352005 1.391222669 2.812491699 -0.079194981 0.793992923
## [246] 0.160707278 0.856195800 1.028701092 0.974676599 0.718697669
## [251] 1.318808130 -0.836785982 1.926158338 0.311743554 -0.877331911
## [256] 1.823077978 -0.605523308 0.484029493 0.451002857 1.216328192
## [261] 0.141117129 0.594160932 -0.300799512 0.218927060 -0.056198180
## [266] -0.498120501 -0.664460833 -0.205382853 1.169024087 0.671026347
## [271] -0.398818264 1.355470314 0.764150938 0.992377247 -0.038016008
## [276] -0.110930387 -1.460641147 -0.208263843 0.207574817 -0.836785982
## [281] -1.235140358 -1.508093443 -0.918081219 -0.759286885 0.757378365
## [286] -0.892969903 -0.424766489 -1.127485479 -0.239381866 -1.513017996
## [291] -0.934132873 -1.047357068 -0.886146303 -1.038550292 -1.460641147
## [296] -1.101432806 -0.253440050 -1.127485479 -1.038550292 -0.909056676
## [301] -1.163171933 -0.436695075 -1.064905886 -0.360329507 -0.556049735
## [306] -1.211192661 -1.627693317 -0.858054530 -0.686613784 -0.571478882
## [311] -0.527921566 -0.715914034 0.746107306 -0.723876456 0.387238121
## [316] -0.723876456 -0.270815112 -1.064905886 -0.694276487 -1.462044259
## [321] -1.101432806 -0.842279246 -0.909056676
## attr(,"gstari")
## [1] FALSE
## attr(,"call")
## localG(x = mpsz_demand_selected_sp$eldercare_count, listw = dnb_lw)
## attr(,"class")
## [1] "localG"
eldercare.gi <- cbind(mpsz_demand_selected_sp, as.matrix(gi.fixed))
names(eldercare.gi)[15] <- "gstat"
#colnames(eldercare.gi)[23] <- "gstat"
tm_shape(eldercare.gi) +
tm_fill(col = "gstat",
style = "pretty",
palette="-RdBu",
title = "local Gi") +
tm_borders(alpha = 0.5)
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## Variable(s) "gstat" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
Kernel Density Maps
CHAS Clinics
bindingbox <- st_bbox(sg)
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
owin <- as(sg, "owin")
spatialpoint <- as(sg, "SpatialPolygons")
chas_ppp <- as(chas_spatialpoint, "ppp")
chas_ppp_jit <- rjitter(chas_ppp, retry=TRUE, nsim=1, drop=TRUE)
ppp_chas <- chas_ppp_jit[owin]
selected_osm <- read_osm(bindingbox, ext=1.1)
msgtext <- "KDE Map"
ppp.km <- rescale(ppp_chas, 1000, "km")
kde <- density(ppp.km, sigma=bw.diggle, edge=TRUE, kernel="gaussian")
## Warning: point-in-polygon test had difficulty with 215 points (total score not 0
## or 1)
gridded_kde_bw <- as.SpatialGridDataFrame.im(kde)
kde_bw_raster <- raster(gridded_kde_bw)
projection(kde_bw_raster) <- crs("+init=EPSG:3414 +datum=WGS84 +units=km")
# Plot kernel density map on openstreetmap
tm_shape(selected_osm)+
tm_rgb()+
tm_shape(kde_bw_raster) +
tm_raster("v", alpha=0.5,
palette = "YlOrRd")
Silver Infocomm Junctions
bindingbox <- st_bbox(sg)
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
owin <- as(sg, "owin")
spatialpoint <- as(sg, "SpatialPolygons")
infocomm_ppp <- as(infocomm_spatialpoint, "ppp")
infocomm_ppp_jit <- rjitter(infocomm_ppp, retry=TRUE, nsim=1, drop=TRUE)
ppp_infocomm <- infocomm_ppp_jit[owin]
selected_osm <- read_osm(bindingbox, ext=1.1)
msgtext <- "KDE Map"
ppp.km <- rescale(ppp_infocomm, 1000, "km")
kde <- density(ppp.km, sigma=bw.diggle, edge=TRUE, kernel="gaussian")
## Warning: point-in-polygon test had difficulty with 215 points (total score not 0
## or 1)
## Warning: point-in-polygon test had difficulty with 215 points (total score not 0
## or 1)
gridded_kde_bw <- as.SpatialGridDataFrame.im(kde)
kde_bw_raster <- raster(gridded_kde_bw)
projection(kde_bw_raster) <- crs("+init=EPSG:3414 +datum=WGS84 +units=km")
# Plot kernel density map on openstreetmap
tm_shape(selected_osm)+
tm_rgb()+
tm_shape(kde_bw_raster) +
tm_raster("v", alpha=0.5,
palette = "YlOrRd")
Eldercare Services
bindingbox <- st_bbox(sg)
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
owin <- as(sg, "owin")
spatialpoint <- as(sg, "SpatialPolygons")
eldercare_ppp <- as(eldercare_spatialpoint, "ppp")
eldercare_ppp_jit <- rjitter(eldercare_ppp, retry=TRUE, nsim=1, drop=TRUE)
ppp_eldercare <- eldercare_ppp_jit[owin]
selected_osm <- read_osm(bindingbox, ext=1.1)
msgtext <- "KDE Map"
ppp.km <- rescale(ppp_eldercare, 1000, "km")
kde <- density(ppp.km, sigma=bw.diggle, edge=TRUE, kernel="gaussian")
## Warning: point-in-polygon test had difficulty with 215 points (total score not 0
## or 1)
## Warning: point-in-polygon test had difficulty with 215 points (total score not 0
## or 1)
gridded_kde_bw <- as.SpatialGridDataFrame.im(kde)
kde_bw_raster <- raster(gridded_kde_bw)
projection(kde_bw_raster) <- crs("+init=EPSG:3414 +datum=WGS84 +units=km")
# Plot kernel density map on openstreetmap
tm_shape(selected_osm)+
tm_rgb()+
tm_shape(kde_bw_raster) +
tm_raster("v", alpha=0.5,
palette = "YlOrRd")
